{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import lu\n",
    "a=np.array([[3,-1,2,4],[1,2,3,4],[2,-2,-1,4],[1,2,3,4]],dtype='float32')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LU分解过程，先计算u的第一行（等于a的第一行），在计算l的第一列（等于a[：，0]/u[0][0]）\n",
    "#再逐列计算，算u（kj）,再算l(ik),再算下一列k+1\n",
    "#计算u(2,2)\n",
    "def myLU(a):\n",
    "    row,col=a.shape\n",
    "    u=np.zeros((row,row),dtype='float32')\n",
    "    l=np.identity((row),dtype='float32')    \n",
    "    u[0]=a[0]\n",
    "    l[:,0]=a[:,0]/u[0][0]\n",
    "    \n",
    "    def calU(k,j):\n",
    "        tempU=0\n",
    "        for i in range(0,k):\n",
    "            tempU=l[k][i]*u[i][j]+tempU\n",
    "        tempU=a[k][j]-tempU\n",
    "        return tempU\n",
    "\n",
    "    def calL(k,j):\n",
    "        tempL=0\n",
    "        for i in range(0,k):\n",
    "            tempL=l[k][i]*u[i][j]+tempL\n",
    "        tempL=(a[k][j]-tempL)/u[j][j]\n",
    "        return tempL\n",
    "\n",
    "    for i in range(1,row):\n",
    "        for j in range(i,row):\n",
    "            u[i][j]=calU(i,j)\n",
    "        for j in range(i+1,row):\n",
    "            l[j][i]=calL(j,i)\n",
    "    return l,u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00]\n",
      " [ 3.3333334e-01  1.0000000e+00  0.0000000e+00  0.0000000e+00]\n",
      " [ 6.6666669e-01 -5.7142860e-01  1.0000000e+00  0.0000000e+00]\n",
      " [ 3.3333334e-01  1.0000000e+00 -5.9604645e-08  1.0000000e+00]]\n",
      "[[ 3.0000000e+00 -1.0000000e+00  2.0000000e+00  4.0000000e+00]\n",
      " [ 0.0000000e+00  2.3333333e+00  2.3333333e+00  2.6666665e+00]\n",
      " [ 0.0000000e+00  0.0000000e+00 -1.0000000e+00  2.8571429e+00]\n",
      " [ 0.0000000e+00  0.0000000e+00  0.0000000e+00  2.8950828e-07]]\n"
     ]
    }
   ],
   "source": [
    "l,u=myLU(a)\n",
    "print(l)\n",
    "print(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.       , -1.       ,  2.       ,  4.       ],\n",
       "       [ 1.       ,  1.9999999,  3.       ,  4.       ],\n",
       "       [ 2.       , -2.       , -1.       ,  4.       ],\n",
       "       [ 1.       ,  1.9999999,  3.       ,  4.       ]], dtype=float32)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l.dot(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[1., 0., 0., 0.],\n",
       "        [0., 1., 0., 0.],\n",
       "        [0., 0., 1., 0.],\n",
       "        [0., 0., 0., 1.]], dtype=float32),\n",
       " array([[ 1.        ,  0.        ,  0.        ,  0.        ],\n",
       "        [ 0.33333334,  1.        ,  0.        ,  0.        ],\n",
       "        [ 0.6666667 , -0.57142854,  1.        ,  0.        ],\n",
       "        [ 0.33333334,  1.        , -0.        ,  1.        ]],\n",
       "       dtype=float32),\n",
       " array([[ 3.       , -1.       ,  2.       ,  4.       ],\n",
       "        [ 0.       ,  2.3333333,  2.3333333,  2.6666665],\n",
       "        [ 0.       ,  0.       , -1.0000001,  2.8571424],\n",
       "        [ 0.       ,  0.       ,  0.       ,  0.       ]], dtype=float32))"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lu(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
